{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4a02d642-5482-447b-a830-1bb571b0c1c8",
   "metadata": {},
   "outputs": [],
   "source": [
    "using LinearAlgebra\n",
    "using Random\n",
    "using Plots\n",
    "using Statistics \n",
    "using DataFrames\n",
    "using DelimitedFiles\n",
    "\n",
    "# 定义参数范围和步长\n",
    "J1_start = 1.0  # J1 的起始值\n",
    "J1_end = 2.0    # J1 的结束值\n",
    "J1_step = 0.5   # J1 的步长\n",
    "J1_list = collect(range(J1_start, stop=J1_end, step=J1_step))\n",
    "\n",
    "\n",
    "J2_list = [-2]\n",
    "\n",
    "\n",
    "J3_list = [0]\n",
    "\n",
    "\n",
    "U_list = [1]\n",
    "\n",
    "\n",
    "mu_list = [0.5]\n",
    "\n",
    "for J1 in J1_list\n",
    "    for J2 in J2_list\n",
    "        for J3 in J3_list\n",
    "            for U in U_list\n",
    "                for mu in mu_list\n",
    "                    n = 4\n",
    "                    N = n * n\n",
    "\n",
    "                    function generate(n)\n",
    "                        return zeros(Int, n*n)\n",
    "                    end\n",
    "\n",
    "                    # Function to compute first, second, and third nearest neighbors\n",
    "                    function neighbor_honeycomb(n)\n",
    "                        nbor1 = zeros(Int, 3, n * n)\n",
    "                        nbor2 = zeros(Int, 6, n * n)\n",
    "                        nbor3 = zeros(Int, 3, n * n)\n",
    "\n",
    "                        for ispin in 1:(n * n)\n",
    "                            iy = div(ispin - 1, n)\n",
    "                            ix = ispin - 1 - iy * n\n",
    "\n",
    "                            # Calculate indices for neighboring spins with periodic boundary conditions\n",
    "                            ixp = mod(ix + 1, n)\n",
    "                            iyp = mod(iy + 1, n)\n",
    "                            ixm = mod(ix - 1 + n, n)\n",
    "                            iym = mod(iy - 1 + n, n)\n",
    "\n",
    "                            # Calculate second nearest indices for neighboring spins\n",
    "                            ixpp = mod(ix + 2, n)\n",
    "                            ixmm = mod(ix - 2 + n, n)\n",
    "                            iypp = mod(iy + 2, n)\n",
    "                            iymm = mod(iy - 2 + n, n)\n",
    "\n",
    "                            # First nearest neighbors\n",
    "                            nbor1[1, ispin] = iy * n + ixp + 1\n",
    "                            nbor1[3, ispin] = iy * n + ixm + 1\n",
    "\n",
    "                            if iy % 2 == 1\n",
    "                                if ix % 2 == 1\n",
    "                                    nbor1[2, ispin] = iyp * n + ix + 1\n",
    "\n",
    "                                    # Second nearest neighbors\n",
    "                                    nbor2[1, ispin] = iy  * n + ixpp + 1\n",
    "                                    nbor2[2, ispin] = iyp * n + ixp  + 1\n",
    "                                    nbor2[3, ispin] = iyp * n + ixm  + 1\n",
    "                                    nbor2[4, ispin] = iy  * n + ixmm + 1\n",
    "                                    nbor2[5, ispin] = iym * n + ixm  + 1\n",
    "                                    nbor2[6, ispin] = iym * n + ixp  + 1\n",
    "\n",
    "                                    # Third nearest neighbors\n",
    "                                    nbor3[1, ispin] = mod(iyp * n + ixpp, n*n) + 1\n",
    "                                    nbor3[2, ispin] = mod(iyp * n + ixmm, n*n) + 1\n",
    "                                    nbor3[3, ispin] = mod(iym * n + ix,   n*n) + 1\n",
    "                                else\n",
    "                                    nbor1[2, ispin] = iym * n + ix + 1\n",
    "\n",
    "                                    # Second nearest neighbors\n",
    "                                    nbor2[1, ispin] = iy  * n + ixpp + 1\n",
    "                                    nbor2[2, ispin] = iyp * n + ixp  + 1\n",
    "                                    nbor2[3, ispin] = iyp * n + ixm  + 1\n",
    "                                    nbor2[4, ispin] = iy  * n + ixmm + 1\n",
    "                                    nbor2[5, ispin] = iym * n + ixm  + 1\n",
    "                                    nbor2[6, ispin] = iym * n + ixp  + 1\n",
    "\n",
    "                                    # Third nearest neighbors\n",
    "                                    nbor3[1, ispin] = mod(iyp * n + ix, n*n) + 1\n",
    "                                    nbor3[2, ispin] = mod(iym * n + ixmm, n*n) + 1\n",
    "                                    nbor3[3, ispin] = mod(iym * n + ixpp, n*n) + 1\n",
    "                                end\n",
    "                            else\n",
    "                                if ix % 2 == 1\n",
    "                                    nbor1[2, ispin] = iym * n + ix + 1\n",
    "\n",
    "                                    # Second nearest neighbors\n",
    "                                    nbor2[1, ispin] = iy  * n + ixpp + 1\n",
    "                                    nbor2[2, ispin] = iyp * n + ixp  + 1\n",
    "                                    nbor2[3, ispin] = iyp * n + ixm  + 1\n",
    "                                    nbor2[4, ispin] = iy  * n + ixmm + 1\n",
    "                                    nbor2[5, ispin] = iym * n + ixm  + 1\n",
    "                                    nbor2[6, ispin] = iym * n + ixp  + 1\n",
    "\n",
    "                                    # Third nearest neighbors\n",
    "                                    nbor3[1, ispin] = mod(iyp * n + ix, n*n) + 1\n",
    "                                    nbor3[2, ispin] = mod(iym * n + ixmm, n*n) + 1\n",
    "                                    nbor3[3, ispin] = mod(iym * n + ixpp, n*n) + 1\n",
    "                                else\n",
    "                                    nbor1[2, ispin] = iyp * n + ix + 1\n",
    "\n",
    "                                    # Second nearest neighbors\n",
    "                                    nbor2[1, ispin] = iy  * n + ixpp + 1\n",
    "                                    nbor2[2, ispin] = iyp * n + ixp  + 1\n",
    "                                    nbor2[3, ispin] = iyp * n + ixm  + 1\n",
    "                                    nbor2[4, ispin] = iy  * n + ixmm + 1\n",
    "                                    nbor2[5, ispin] = iym * n + ixm  + 1\n",
    "                                    nbor2[6, ispin] = iym * n + ixp  + 1\n",
    "\n",
    "                                    # Third nearest neighbors\n",
    "                                    nbor3[1, ispin] = mod(iyp * n + ixpp, n*n) + 1\n",
    "                                    nbor3[2, ispin] = mod(iyp * n + ixmm, n*n) + 1\n",
    "                                    nbor3[3, ispin] = mod(iym * n + ix, n*n) + 1\n",
    "                                end\n",
    "                            end\n",
    "                        end\n",
    "                        return nbor1, nbor2, nbor3\n",
    "                    end\n",
    "\n",
    "                    function calculate(ising, nbor1, nbor2, nbor3, J1, J2, J3, N, U, mu)\n",
    "                        e = 0.0\n",
    "                        neighbor_sum = 0.0\n",
    "                        for i in 1:N\n",
    "                            for j in 1:3\n",
    "                                if ising[i] != ising[nbor1[j,i]]\n",
    "                                    neighbor_sum = neighbor_sum - J1\n",
    "                                end\n",
    "                            end\n",
    "                            for j in 1:6\n",
    "                                if ising[i] != ising[nbor2[j, i]]\n",
    "                                    neighbor_sum = neighbor_sum - J2\n",
    "                                end\n",
    "                            end\n",
    "                            for j in 1:3\n",
    "                                if ising[i] != ising[nbor3[j, i]]\n",
    "                                    neighbor_sum = neighbor_sum - J3\n",
    "                                end\n",
    "                            end\n",
    "                        end\n",
    "                        e = e + neighbor_sum\n",
    "                        for i in 1:N\n",
    "                            e = e + (U / 2) * ising[i] * (ising[i] - 1)\n",
    "                        end\n",
    "                        for i in 1:N\n",
    "                            e = e - mu * ising[i]\n",
    "                        end\n",
    "                        e = e / 2\n",
    "                        m = (sum(ising))\n",
    "                        return e, m\n",
    "                    end\n",
    "\n",
    "                    ising = generate(n)\n",
    "                    nbor1, nbor2, nbor3 = neighbor_honeycomb(n)\n",
    "\n",
    "                    e_steps = zeros(3^N)\n",
    "                    m_steps = zeros(3^N)\n",
    "                    ising_steps=zeros(3^N, N)\n",
    "\n",
    "                    for k1 in 0:3^N - 1\n",
    "                        num1 = lpad(string(k1, base=3), N, '0')\n",
    "                        for j in 1:N\n",
    "                            ising[j] = parse(Int, num1[j])\n",
    "                        end    \n",
    "                        ising_steps[k1+1,:]=ising[:]\n",
    "                        e_steps[k1 + 1], m_steps[k1 + 1] = calculate(ising, nbor1, nbor2, nbor3, J1, J2, J3, N, U, mu)\n",
    "                    end\n",
    "\n",
    "                    min_energy, index = findmin(e_steps)\n",
    "                    println(\"最小能量 (min energy): \", min_energy)\n",
    "\n",
    "                    ising_steps[index,:]\n",
    "                    configuration = ising_steps[index, :]\n",
    "\n",
    "                    count_min_energy = count(x -> x == min_energy, e_steps)\n",
    "                    indices = findall(x -> x == min_energy, e_steps)\n",
    "                    min_configurations = ising_steps[indices, :] \n",
    "\n",
    "                    # 给定的坐标信息\n",
    "                    x = [0.0, 0.8660254037844386, 1.7320508075688772, 2.598076211353316, 0.0, 0.8660254037844386, 1.7320508075688772, 2.598076211353316, 0.0, 0.8660254037844386, 1.7320508075688772, 2.598076211353316, 0.0, 0.8660254037844386, 1.7320508075688772, 2.598076211353316]\n",
    "                    y = [0.0, -0.5, 0.0, -0.5, 1.0, 1.5, 1.0, 1.5, 3.0, 2.5, 3.0, 2.5, 4.0, 4.5, 4.0, 4.5]\n",
    "\n",
    "                    # 颜色映射：0 映射为蓝色，1 映射为红色，2 映射为绿色\n",
    "                    color_mapping = Dict(0.0 => :blue, 1.0 => :red, 2.0 => :green)\n",
    "\n",
    "                    folder_name = \"J1_$(J1)_J2_$(J2)_J3_$(J3)_U_$(U)_μ_$(mu)\"\n",
    "                    ispath(folder_name) || mkpath(folder_name)\n",
    "\n",
    "                    # 循环绘制每个与最小能量对应的自旋配置的散点图\n",
    "                    for i in 1:count_min_energy \n",
    "                        p = plot()\n",
    "                        z = min_configurations[i, :]\n",
    "                        colors = [color_mapping[z_i] for z_i in z]\n",
    "                        scatter!(x, y, color = colors, aspect_ratio = 1)\n",
    "                        display(p)\n",
    "                        savefig(p, joinpath(folder_name, \"simple_plot_$i.png\"))\n",
    "                    end\n",
    "                end\n",
    "            end\n",
    "        end\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6b9c7e41-d562-4273-bb03-9bd56110c3ce",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.11.4",
   "language": "julia",
   "name": "julia-1.11"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
